Summary

We obtain the IRFs for some of the best models (with respect to AIC/BIC/Shapiro-Wilk/Jarque-Bera/Ljung-Box). We use either estimated values for the static shock transmission matrix \(B\) or impose the long-run restriction suggested by Blanchard and Quah (1989).

In addition, we generate the IRFs of Blanchard and Quah (1989) and GMR.

Preliminaries

Load and plot data

tt = readRDS(paste0(params$PATH, params$JOBID, "/",
                    "e_residualcheck_end.rds")) %>% arrange(rk_bic)

DATASET = readRDS(params$PATH_DATASET)
dygraphs::dygraph(DATASET)
DATASET = DATASET %>% as.matrix()
DIM_OUT = dim(DATASET)[2]
N_OBS = dim(DATASET)[1]

Define some functions for IRFs

# Permute static shock transmission matrix B and change signs

permute_chgsign = function(irf_array, 
                           perm = rep(1, dim(irf_array)[1]), 
                           sign = rep(1, dim(irf_array)[1])){
  
  dim_out = dim(irf_array)[1]
  
  perm_mat = diag(dim_out)
  perm_mat = perm_mat[,perm]
  
  sign = diag(sign)
  
  ll = map(1:dim(irf_array)[3], ~ irf_array[,,.x] %*% perm_mat %*% sign)
  
  irf_array = array(0, dim = c(dim_out, dim_out, dim(irf_array)[3]))
  for (ix in 1:dim(irf_array)[3]){
    irf_array[,,ix] = ll[[ix]] 
  }
  
  return(irf_array)
}

# Plot IRF

plot_irf = function(irf_array){
  plot(pseries(irf_array %>% polm(), lag.max = 40))
}


# Calculate the cumulative sum of the IRF coefficients for a given element 

irf_cumsum = function(irf_array, el_mat){
  
  stopifnot("*el_mat* should be a matrix with two columns where each row contains the indices for which the cumulative sum should be calculated" = 
            dim(el_mat)[2] == 2)
  
  for (ix in (1:nrow(el_mat))){
    irf_array[el_mat[ix,1], el_mat[ix,2], ] = irf_array[el_mat[ix,1], el_mat[ix,2], ] %>% cumsum()
  }
  
  return(irf_array)
}


# Rotate the static shock transmission matrix such that there is no long-run effect

rotate_longrun = function(irf_array){

  dim_out = dim(irf_array)[1]
  
  stopifnot("Number of outputs must be equal to 2" = dim_out == 2)
  
  k1 = freqresp(armamod(lmfd(polm(diag(dim_out)), polm(irf_array)), sigma_L = diag(dim_out)), n.f = 1)$frr
  lq = lq_decomposition(matrix(c(Re(k1)), nrow = dim_out))
  Q = lq$q

  jj = map(1:dim(irf_array)[3], ~ irf_array[,,.x] %*% t(Q))
  
  irf_array = array(0, dim = c(dim_out, dim_out,dim(irf_array)[3]))
  for (ix in 1:dim(irf_array)[3]){
    irf_array[,,ix] = jj[[ix]] 
  }

return(irf_array)
}

Model statistics

tt %>%
  select(-nr, -p_plus_q, -rk_mle, -contains("value"), -contains("cov"), -contains("pval")) %>% 
  filter(normality_flag == 0, lb_flag == 0) %>% 
  dtable()

Blanchard Quah AR(8)

Function that calculates the cumulative IRF for some variables of choice. Here, we are interested in log(GNP) and not its change.

Estimate a VAR(8) model with the svars package.

vars_bq = VAR(DATASET, p = 8, type = "none")
vars_irf_bq = irf(vars_bq, n.ahead = 100)

Transform it to rationalmatrices/RLDM style and create a first plot.

irf_bq = array(0, dim = c(2,2,101))
irf_bq[,1,] = t(vars_irf_bq$irf$rGDPgrowth_demeaned)
irf_bq[,2,] = t(vars_irf_bq$irf$unemp_detrended)

plot(pseries(polm(irf_bq), lag.max = 100))

Calculate the long-run response, LQ decompose it, rotate the static shock impact matrix, transform the transfer function s.t. the demand shock is transitory.

k1 = matrix(c(colSums(vars_irf_bq$irf$rGDPgrowth_demeaned),
              colSums(vars_irf_bq$irf$unemp_detrended)),
            nrow = 2)
lq = lq_decomposition(k1)
Q = lq$q
 
ii = map(1:101, ~ irf_bq[,,.x] %*% t(Q))

irf_bq_lr = array(0, dim = c(2,2,101))
for (ix in 1:101){
  irf_bq_lr[,,ix] = ii[[ix]] 
}
rm(k1, Q, ii, irf_bq)

Check the result : transfer function evaluated at \(1\) should be triangular.

apply(irf_bq_lr, c(1,2), sum) %>% zapsmall()
##           [,1]     [,2]
## [1,]  0.622691 0.000000
## [2,] -0.161077 5.558043

Plot the rotated IRF

plot(pseries(irf_bq_lr %>% polm(), lag.max = 40))

irf_bq_lr_cumsum = irf_cumsum(irf_bq_lr, matrix(c(1,1,1,2), nrow = 2))

Plot the rotated IRF where the response to (log) GNP is cumulated.

plot(pseries(irf_bq_lr_cumsum %>% polm(), lag.max = 40))

Permute columns and change sign such that the plots align with the ones in Blanchard and Quah.

jj = map(1:101, ~ irf_bq_lr_cumsum[,,.x] %*% matrix(c(0,-1,1,0), nrow = 2))

irf_bq_lr_cumsum_permuted = array(0, dim = c(2,2,101))
for (ix in 1:101){
  irf_bq_lr_cumsum_permuted[,,ix] = jj[[ix]] 
}

plot(pseries(irf_bq_lr_cumsum_permuted %>% polm(), lag.max = 40))

GMR VARMA(4,1, k = 1)

if(!file.exists("../local_data/gmr_results_bq.rds")){
  load("~/r_projects/gmr_ssvarma/results/BQ/results.MLE.4lags.res")
  write_rds(res.estim.MLE, "../local_data/gmr_results_bq.rds")
  gmr_results_bq = res.estim.MLE
} else{
  gmr_results_bq = read_rds("../local_data/gmr_results_bq.rds")
}

Convert GMR results to rationalmatrices/RLDM style.

ar_gmr = gmr_results_bq$Phi
ma_gmr = gmr_results_bq$Theta
B_gmr = gmr_results_bq$C

polm_ar = array(c(c(diag(2)),
                  -c(ar_gmr)),
                dim = c(2,2,5)) %>% polm()

polm_ma = array(c(c(diag(2)),
                  -c(ma_gmr)),
                dim = c(2,2,2)) %>% polm()

B = gmr_results_bq$C

lmfd_obj = lmfd(polm_ar, polm_ma %r% B) 

Check poles and zeros

polm_ar %>% zeroes() %>% abs()
## [1] 1.192445 1.192445 1.324240 1.447863 1.447863 2.216027 2.216027 4.102746
polm_ma %>% zeroes() %>% abs()
## [1] 0.565828 2.519183
irf_gmr = pseries(lmfd_obj, lag.max = 40) %>% unclass()
plot_irf(irf_gmr)

irf_gmr_cumsum = irf_cumsum(irf_gmr, 
                           el_mat = matrix(c(1,1,1,2), nrow = 2))

plot_irf(irf_gmr_cumsum)

Permute columns such that results align

irf_gmr_cumsum_permuted = permute_chgsign(irf_gmr_cumsum, perm = c(2,1))
plot_irf(irf_gmr_cumsum_permuted)

BF VARMA(1,2, kappa = 1, k = 0)

Extract appropriate inter-valued parameters.

tt %>% 
  filter(p==1 & q==2 & kappa==1 & k==0) %>% 
  select(p_plus_q, p,q, kappa, k, starts_with("rk_"), cov_el_sum, value_final, value_aic, value_bic, everything()) -> tt_bf_12_2 

Create an ARMA-WHF model:

bf_params_12_2 = tt_bf_12_2$params_deep_final %>% .[[1]]
bf_tmpl_12_2 = tt_bf_12_2$tmpl %>% .[[1]]

write_rds(bf_params_12_2, "../local_data/p_whf/best_mod_rob/bf_params_12_2.rds")
write_rds(bf_tmpl_12_2, "../local_data/p_whf/best_mod_rob/bf_tmpl_12_2.rds")


bf_armawhfmod_12_2 = fill_tmpl_whf_rev(bf_params_12_2, bf_tmpl_12_2)

Create IRF, cumulate GNP, and plot:

bf_irf_12_2 = irf_whf(bf_params_12_2, bf_tmpl_12_2, n_lags = 40) %>% unclass()
bf_irf_12_2_cumsum = bf_irf_12_2 %>% irf_cumsum(el_mat = matrix(c(1,1,1,2), nrow = 2))
bf_irf_12_2_cumsum %>% plot_irf()

bf_irf_12_2_cumsum_permuted = bf_irf_12_2_cumsum

Rotate such that BQ long-run restriction is satisfied, plot, sign-permute, and plot again.

bf_irf_12_2_lr = bf_irf_12_2 %>% rotate_longrun()
bf_irf_12_2_lr_cumsum = bf_irf_12_2_lr %>% irf_cumsum(el_mat = matrix(c(1,1,1,2), nrow = 2))
bf_irf_12_2_lr_cumsum %>% plot_irf()

bf_irf_12_2_lr_cumsum_permuted = bf_irf_12_2_lr_cumsum %>% permute_chgsign(perm = c(2,1),
                                                                           sign = c(-1,1))
bf_irf_12_2_lr_cumsum_permuted %>% plot_irf()

BF VARMA(1,3, kappa = 1, k = 0)

Extract appropriate inter-valued parameters.

tt %>% 
  filter(p==1 & q==3 & kappa==1 & k==0) %>% 
  select(p_plus_q, p,q, kappa, k, starts_with("rk_"), cov_el_sum, value_final, value_aic, value_bic, everything()) -> tt_bf_13_2 

Create an ARMA-WHF model:

bf_params_13_2 = tt_bf_13_2$params_deep_final %>% .[[1]]
bf_tmpl_13_2 = tt_bf_13_2$tmpl %>% .[[1]]

bf_armawhfmod_13_2 = fill_tmpl_whf_rev(bf_params_13_2, bf_tmpl_13_2)

Create IRF, cumulate GNP, and plot:

bf_irf_13_2 = irf_whf(bf_params_13_2, bf_tmpl_13_2, n_lags = 40) %>% unclass()
bf_irf_13_2_cumsum = bf_irf_13_2 %>% irf_cumsum(el_mat = matrix(c(1,1,1,2), nrow = 2))
bf_irf_13_2_cumsum %>% plot_irf()

bf_irf_13_2_cumsum_permuted = bf_irf_13_2_cumsum

Rotate such that BQ long-run restriction is satisfied, plot, sign-permute, and plot again.

bf_irf_13_2_lr = bf_irf_13_2 %>% rotate_longrun()
bf_irf_13_2_lr_cumsum = bf_irf_13_2_lr %>% irf_cumsum(el_mat = matrix(c(1,1,1,2), nrow = 2))
bf_irf_13_2_lr_cumsum %>% plot_irf()

bf_irf_13_2_lr_cumsum_permuted = bf_irf_13_2_lr_cumsum %>% permute_chgsign(perm = c(2,1),
                                                                           sign = c(-1,1))
bf_irf_13_2_lr_cumsum_permuted %>% plot_irf()

BF VARMA(3,1, kappa = 1, k = 0)

Extract appropriate inter-valued parameters.

tt %>% 
  filter(p==3 & q==1 & kappa==1 & k==0) %>% 
  select(p_plus_q, p,q, kappa, k, starts_with("rk_"), cov_el_sum, value_final, value_aic, value_bic, everything()) -> tt_bf_31_2 

Create an ARMA-WHF model:

bf_params_31_2 = tt_bf_31_2$params_deep_final %>% .[[1]]
bf_tmpl_31_2 = tt_bf_31_2$tmpl %>% .[[1]]

bf_armawhfmod_31_2 = fill_tmpl_whf_rev(bf_params_31_2, bf_tmpl_31_2)

Create IRF, cumulate GNP, and plot:

bf_irf_31_2 = irf_whf(bf_params_31_2, bf_tmpl_31_2, n_lags = 40) %>% unclass()
bf_irf_31_2_cumsum = bf_irf_31_2 %>% irf_cumsum(el_mat = matrix(c(1,1,1,2), nrow = 2))
bf_irf_31_2_cumsum %>% plot_irf()

Sign-permute such that shocks are labelled consistently.

bf_irf_31_2_cumsum_permuted = bf_irf_31_2_cumsum %>% permute_chgsign(perm = c(2,1),
                                                                     sign = c(1,1))
bf_irf_31_2_cumsum_permuted %>% plot_irf()

Rotate such that BQ long-run restriction is satisfied, plot, sign-permute, and plot again.

bf_irf_31_2_lr = bf_irf_31_2 %>% rotate_longrun()
bf_irf_31_2_lr_cumsum = bf_irf_31_2_lr %>% irf_cumsum(el_mat = matrix(c(1,1,1,2), nrow = 2))
bf_irf_31_2_lr_cumsum %>% plot_irf()

bf_irf_31_2_lr_cumsum_permuted = bf_irf_31_2_lr_cumsum %>% permute_chgsign(perm = c(2,1),
                                                                           sign = c(-1,1))
bf_irf_31_2_lr_cumsum_permuted %>% plot_irf()

BF VARMA(1,1, kappa = 1)

Extract appropriate inter-valued parameters.

tt %>% 
  filter(p==1 & q==1 & kappa==1 & k==0) %>% 
  select(p_plus_q, p,q, kappa, k, starts_with("rk_"), cov_el_sum, value_final, value_aic, value_bic, everything()) -> tt_bf_11_2 

Create an ARMA-WHF model:

bf_params_11_2 = tt_bf_11_2$params_deep_final %>% .[[1]]
bf_tmpl_11_2 = tt_bf_11_2$tmpl %>% .[[1]]

bf_armawhfmod_11_2 = fill_tmpl_whf_rev(bf_params_11_2, bf_tmpl_11_2)

Create IRF, cumulate GNP, and plot:

bf_irf_11_2 = irf_whf(bf_params_11_2, bf_tmpl_11_2, n_lags = 40) %>% unclass()
bf_irf_11_2_cumsum = bf_irf_11_2 %>% irf_cumsum(el_mat = matrix(c(1,1,1,2), nrow = 2))
bf_irf_11_2_cumsum %>% plot_irf()

Sign-permute such that shocks are labelled consistently.

bf_irf_11_2_cumsum_permuted = bf_irf_11_2_cumsum %>% permute_chgsign(perm = c(2,1),
                                                                     sign = c(-1,-1))
bf_irf_11_2_cumsum_permuted %>% plot_irf()

Rotate such that BQ long-run restriction is satisfied, plot, sign-permute, and plot again.

bf_irf_11_2_lr = bf_irf_11_2 %>% rotate_longrun()
bf_irf_11_2_lr_cumsum = bf_irf_11_2_lr %>% irf_cumsum(el_mat = matrix(c(1,1,1,2), nrow = 2))
bf_irf_11_2_lr_cumsum %>% plot_irf()

bf_irf_11_2_lr_cumsum_permuted = bf_irf_11_2_lr_cumsum %>% permute_chgsign(perm = c(2,1),
                                                                           sign = c(-1,1))
bf_irf_11_2_lr_cumsum_permuted %>% plot_irf()

BF VARMA(4,1, kappa = 0, k = 1): GMR’s model

Extract appropriate inter-valued parameters.

tt %>% 
  filter(p==4 & q==1 & kappa==0 & k==1) %>% 
  select(p_plus_q, p,q, kappa, k, starts_with("rk_"), cov_el_sum, value_final, value_aic, value_bic, everything()) -> tt_bf_41_1 

Create an ARMA-WHF model:

bf_params_41_1 = tt_bf_41_1$params_deep_final %>% .[[1]]
bf_tmpl_41_1 = tt_bf_41_1$tmpl %>% .[[1]]

bf_armawhfmod_41_1 = fill_tmpl_whf_rev(bf_params_41_1, bf_tmpl_41_1)

Create IRF, cumulate GNP, and plot:

bf_irf_41_1 = irf_whf(bf_params_41_1, bf_tmpl_41_1, n_lags = 40) %>% unclass()
bf_irf_41_1_cumsum = bf_irf_41_1 %>% irf_cumsum(el_mat = matrix(c(1,1,1,2), nrow = 2))
bf_irf_41_1_cumsum %>% plot_irf()

Sign-permute such that shocks are labelled consistently.

bf_irf_41_1_cumsum_permuted = bf_irf_41_1_cumsum %>% permute_chgsign(perm = c(2,1))
bf_irf_41_1_cumsum_permuted %>% plot_irf()

Rotate such that BQ long-run restriction is satisfied, plot, sign-permute, and plot again.

bf_irf_41_1_lr = bf_irf_41_1 %>% rotate_longrun()
bf_irf_41_1_lr_cumsum = bf_irf_41_1_lr %>% irf_cumsum(el_mat = matrix(c(1,1,1,2), nrow = 2))
bf_irf_41_1_lr_cumsum %>% plot_irf()

bf_irf_41_1_lr_cumsum_permuted = bf_irf_41_1_lr_cumsum %>% permute_chgsign(perm = c(2,1),
                                                                           sign = c(-1,1))
bf_irf_41_1_lr_cumsum_permuted %>% plot_irf()

Summary

list_mods = tibble(BQ_AR8 = irf_bq_lr_cumsum_permuted[,,1:41] %>% list(),
                   GMR_ARMA41_1 = irf_gmr_cumsum_permuted %>% list(),
                   BF_ARMA11_2 = bf_irf_11_2_cumsum_permuted %>% list(),
                   BF_ARMA11_2_lr = bf_irf_11_2_lr_cumsum_permuted %>% list(),
                   BF_ARMA12_2 = bf_irf_12_2_cumsum_permuted %>% list(),
                   BF_ARMA12_2_lr = bf_irf_12_2_lr_cumsum_permuted %>% list(),
                   BF_ARMA13_2 = bf_irf_13_2_cumsum_permuted %>% list(),
                   BF_ARMA13_2_lr = bf_irf_13_2_lr_cumsum_permuted %>% list(),
                   BF_ARMA31_2 = bf_irf_31_2_cumsum_permuted %>% list(),
                   BF_ARMA31_2_lr = bf_irf_31_2_lr_cumsum_permuted %>% list(),
                   BF_ARMA41_1 = bf_irf_41_1_cumsum_permuted %>% list(),
                   BF_ARMA41_1_lr = bf_irf_41_1_lr_cumsum_permuted %>% list())

list_mods %>% pivot_longer(everything()) %>% 
  mutate(Demand2GNP = map(value, ~.x[1,1,]),
         Demand2Unempl = map(value, ~.x[2,1,]),
         Supply2GNP = map(value, ~.x[1,2,]),
         Supply2Unempl = map(value, ~.x[2,2,])) %>% 
  select(-value) %>% 
  pivot_longer(c("Demand2GNP", "Demand2Unempl", "Supply2GNP", "Supply2Unempl"), 
               names_to = "Response_Type", values_to = "Impact") %>% 
  separate(Response_Type, into = c("Shock", "Output"), sep = "2") %>% 
  unnest_longer(Impact, indices_include = TRUE) %>% 
  rename(Lag = Impact_id,
         Model = name) %>% 
  mutate(Lag = Lag-1) -> tibble_irf
tibble_irf %>% 
  # filter(Model != "BF_ARMA11_2") %>% 
  ggplot(aes(x = Lag, y = Impact, color = Model)) +
  geom_line() + geom_point() + facet_grid(Output~Shock) -> ply

plotly::ggplotly(ply)
knitr::knit_exit()